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Abstract 

A new deterministic, numerical method to solve fermion field theories is presented. 
This approach is based on finding solutions Z[J] to the lattice functional equations for 
field theories in the presence of an external source J. Using Grassmann polynomial 
expansions for the generating functional Z, we calculate propagators for systems of in- 
teracting fermions. These calculations are straightforward to perform and are executed 
rapidly compared to Monte Carlo. The bulk of the computation involves a single matrix 
inversion. Because it is not based on a statistical technique, it does not have many of 
the difficulties often encountered when simulating fermions. Since no determinant is 
ever calculated, solutions to problems with dynamical fermions are handled more eas- 
ily. This approach is very flexible, and can be taylored to specific problems based on 
convenience and computational constraints. We present simple examples to illustrate 
the method; more general schemes are desirable for more complicated systems. 
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1 Introduction 



The need for new calculational probes of quantum field theory is well-known. Despite 
the success of various iterative methods including perturbation theory, many problems, 
especially for strongly coupled systems, persist in both high energy and condensed matter 
physics. The advent of high performance computers has made numerical attacks on these 
problems practical. Most of these approaches have roots in Monte Carlo methods. They 
form the basis of many calculations in statistical mechanics, including nonrelativistic field 
theory, and have been co-opted by high energy physics for lattice gauge theory calculations 

0- 

In high energy physics, Monte Carlo methods had remarkable initial successes in cal- 
culating various quantities associated with quantum field theory. Unfortunately, even for 
simple systems, calculations involving dynamical fermions, or performed on large lattices, 
become expensive in both time and computer resources. Despite refined algorithms and 
enormous increases in computer power, progress with Monte Carlo based methods has been 
slow compared to the initial impressive successes. For lattice QCD Q, detailed studies of 
systems with dynamical quarks with low masses may require computer time well beyond 
current capabilities. 

In condensed matter physics, understanding models of strongly correlated electrons us- 
ing numerical methods has serious limitations. Much of this work has been based on exact 
diagonalization or quantum Monte Carlo. Diagonalization is powerful especially for calcu- 
lating dynamical properties, but in practice, is limited to small systems. Typically, 4x4 
clusters are the largest that can be studied. Alternatively, Monte Carlo methods can ex- 
amine larger lattices, but for fermionic systems, simulations face serious problems. This 
includes the so-called "minus sign" problem which makes defining a positive-definite prob- 
ability measure problematic [11]. This difficulty plagues many models, including many of 
those relevant to high T c superconductivity. 

We have formulated a numerical scheme for lattice quantum field theory which does 
not rely on any statistical method and is not restricted to small systems || |], ||, ||. This 
approach is based on considering QFT in the presence of an external source. Using either 
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the operator equations of motion, or equivalently the path integral representation, equations 
can be derived for the vacuum persistence amplitude Z . On a lattice, the result is a set of 
coupled linear differential equations for Z in the discretized sources. For fermionic problems, 
these sources will be anticommuting Grassmann variables. Once we have obtained the 
lattice Green's functions can be extracted by functional differentiation. 

Despite the large number of coupled differential equations for Z, one for each source 
variable, there is nothing unusual about this problem mathematically. We can solve them as 
we do any set of linear differential equations. Since at the end of any calculation the sources 
are set to zero, the obvious first choice is to expand Z in a power series in the sources about 
the origin. For fermionic systems with only a few sites, exact solutions for Z are possible. 
This is due to the self-terminating nature of Grassmann polynomials. In general, we cannot 
solve the theory exactly, and therefore, we must introduce an approximation. For the 
purposes of this paper, we truncate the series expansion for Z after some maximum power 
in the source variables. In this way, we obtain a closed, finite set of linear equations for the 
expansion coefficients accurate to within the limits of the approximation. 

The coupled differential equations have a high degree of symmetry due to the translation 
and reflection invariance of the underlying periodic lattice. This structure can be used to 
great advantage in constructing solutions, greatly reducing the number of equations to be 
considered. This reduced set of algebraic equations for the coefficients of Z is necessarily 
inconsistent; this being an artifact of the approximation. There are many ways to deal with 
this problem. One particularly elegant method is to require that weighted averages of the 
truncated equations vanish. A clever choice of weight functions speeds the convergence of 
the truncated equations to the exact solution. Numerical techniques of this sort are called 
Galerkin methods and are well-known [||]. 

Since our expansion and weight functions depend on the external sources, we call this 
method "the Source Galerkin method" It is different in concept and application from 
any other technique that we know of for solving quantum field theory. As a practical 
matter, fermions and bosons are treated in a similar manner. Consequently, solutions for 
dynamical fermion problems may be possible with existing computers. On the other hand, 
this method is not without difficulties. Principally, the power series grows quickly, making 
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examination of large systems difficult. In these cases, different basis functions must be used. 
The examples presented here illustrate one of the simpler formulations, but not generally 
the most useful. 

This paper is organized as follows. Section Two gives a general overview of the bosonic 
formulation of the Source Galerkin method. Section Three outlines the Source Galerkin 
formulation for fermions. The functional formulation is identical to the boson case except 
that the sources become anticommuting Grassmann variables. The partition function Z 
satisfies a set of coupled Grassmann differential equations. As with bosons, we solve these 
equations by power series in the source variables. In Section Four, we examine small lattice 
systems where exact solutions are possible. For more realistic systems, the series must be 
truncated and a Galerkin procedure devised. We propose a Galerkin method for fermions 
in Section Five and perform calculations on ID self-interacting fermion systems. We close 
by illustrating the flexibility of the Source Galerkin method with respect to the choice of 
expansion functions. We study 2D lattice Gross-Neveu model by examining the leading 
behavior of Z as J — > 0. The first term is equivalent to a Gaussian. We calculate the 
behavior of the chiral condensate as an order parameter for chiral symmetry breaking. A 
method to calculate corrections is outlined. 

2 Source Galerkin Method 

In this section, we review the bosonic formulation of the Source Galerkin method. We 
study quantum field theory in the presence of an external source. The vacuum persistence 
amplitude Z[J\ =j (0j0)j is the generating functional for the Green's functions. It is 
constrained by a functional differential equation. For a self-interacting scalar field theory, 
the dynamics are described by 



Lj(<f>) = 1/2 (a + M 2 )4> + s/4 4 + J$. 



(1) 



Taking vacuum expectation values of the operator equations of motion 



(□+M 2 )(0)j + g{<jfl)j = Jj(0\0)j 



(2) 
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and identifying Euclidean expectation values of fields with functional derivatives of Z[J\ 



f) n 7\ fl 
5J(xi) . . . 5J(x n ) 



yields the functional relation 



, SZ 8 3 Z 



{a + M) JW) +,j jiw = JZ - (4) 

After solving this equation for Z, we can extract all information about our field theory by 
functional differentiation. On a D-dimensional Euclidean lattice, the functional equation 
becomes a set of coupled differential equations 

where the sum is over nearest neighbors. There is one equation per site. 

For finite lattices with N sites, Z as a multivariate function of the N source variables 
J(i). To construct a solution to the differential equations, we can expand Z on any complete 
set of functions in the source variables 

Z = ««0n({J}) (6) 

n 

where a n are the unknowns of the problem. A particularly simple choice is polynomial 
functions. It should be emphasized that other choices are preferred for more complicated 
systems. We consider this possibility in Section Seven. 

The number of independent unknown coefficients can be reduced by exploiting the sym- 
metries of the lattice. For polynomial expansions, Z can be constructed to be invariant 
under the symmetry group of the lattice 

Z = ^ ' flri,m Pn,m (7) 
n,m 

where P n ,m are invariant polynomials in the source variables of order n and with m invari- 
ant classes for a given order. The coefficients a n ^ m are to be determined. The number of 
invariant classes for a given order depends on both the number of lattice sites and on the 
number of symmetry operations for a given lattice. For example, higher dimensional lat- 
tices have larger symmetry groups and therefore have fewer independent unknowns. These 
polynomials form a lattice invariant basis upon which we can construct solutions for Z. 



In order to solve the differential equations, we must specify boundary conditions. Nor- 
malization of the vacuum amplitude implies 



Z[J = 0] = 1 (8) 



and we exclude symmetry breaking 

dZ 



U=o = m)) = 0. (9) 



dJ{i) 

We cannot specify the second derivatives of Z because they correspond to the two-point 
functions. Instead, w truncate Z at some finite order M. The truncation order is dictated 
by computational constraints. As a boundary condition, truncation guarantees that as 
g —* in the interaction theory, the solution reduces to the free field. This is possible only 
for path integral formulations constrained to real integration contours || 0, ]|| . Truncation 
in this sense not only sets a boundary condition, but also introduces an approximation 
scheme. This scheme is made systematic by truncating at successively higher orders, and 
then taking the limit as M goes to infinity. 

We fit our trial solution to the differential equations using the Galerkin method. The 
truncated polynomial is an approximate solution to the differential equations where the 
error due to truncation is called the residual R ||. In order to minimize R, we define an 
inner product in the source space 

(9, f) = £ 9(J(1) ■ ■ ■ J(N))f(J(l) . . . J(N)) [dJ] (10) 

where the integration is over all J(i), and e is considered small. Requiring the inner product 
of R with linearly independent test functions Tk to vanish 

(12, Ti) = 

: (11) 
(R,Tj) = 

generates a set of linear algebraic equations. We can generate as many of these equations 
as we like as long as the Tk are linearly independent. 

In analogy with a variational principle, we chose the test functions as 

n = dJ(T) Pn ' m (12) 



where the P n , m are the invariant polynomial that formed the basis for Zt ■ This choice 
guarantees that we always have as many equations as unknowns, and experience has shown 
that it gives rapid convergence of Zt to the exact solution. We construct as many equations 
as there are independent unknowns in our power series solution These equations can be 
solved by a single matrix inversion where the resulting coefficients are the lattice Green's 
functions. 



3 Fermion Formulation 

We now construct the lattice functional equations for fermions. As with the bosonic case, we 
appeal to a generating function Z to calculate the consequences of our fermionic field theory. 
We begin by considering a ID system of spinless free fermions in the presence of external 
sources. We will see that the inclusion of interactions requires only minor modifications. 
Our system has dynamics described by the following action 

N 1 

s = Yl 2^*) + - w - ^ + M ^(*) + [mm + mm ]. (13) 

i 

where the sum runs over an iV site lattice and periodic boundary conditions are assumed. 
Grassmann sources fj(i),r)(i) live on each site of the lattice and are coupled to the fields. 
Taking the variation of the action with respect to the fields, we obtain the EOM for the 
field operators if>(i) 

+ 1) - ij>(i - 1)] + AT0(i) + r}(i) = 0. (14) 

There are analogous equations for the conjugates. Taking vacuum matrix elements of this 

equation, and identifying expectation values of fields with Grassmann Source derivatives of 

the generating functional Z[fj,rj\. 

f)7 r)7 
' {4>{n)) __ = -<y,(n)), (15) 



drj(i) ' drj(n 

gives the lattice equations for Z. 



SZ - + M^-„ W Z = 0. (16) 



2 dfj(i + 1) dfj(i — 1) dfj(i) 
As with the boson case, there is one equation (and its conjugate) for each site of the lattice. 
Our fermion field theory is completely specified by this set of coupled Grassmann PDEs 



for Z plus boundary conditions. The Green's functions can be calculated as derivatives of 
Z[fj, rj] about 77, 77 = 0. 

We solve these equations with a Grassmann power series 

Z= £ an 1 , ni ...n N ,n N m ni ---v(N) nN va) ni ---v(N) nN . (17) 

n,i,ni=0,l 

Because of Fermi statistics, Z is at most linear in each source variable. This means that 
for finite systems, the power series terminates at a finite order independent of any approx- 
imation such as truncation. Therefore, exact solutions exist for finite series. To reduce the 
number of independent coefficients in this series, we take advantage of symmetries of the 
system. In addition to the space-time symmetries, the Action conserves fermion number 

and is invariant under charge conjugation 

For bipartite lattices, we have the additional symmetries 

— > —ip(i) if i = odd 
■ifj(i) —> if i = even. 

We construct the terms in our series explicitly to reflect all these symmetries. 

4 Exact Solutions for Four Site Model 

To be specific, consider a free fermion field on four sites. We place complex spinless fields 
and sources on each site. This gives eight coupled Grassmann differential equations for Z 
including the conjugate equations. Since they are all related by symmetry operations, it 
is sufficient to use only one equation to constrain our solution as long as Z is constructed 
explicitly to reflects these symmetries. We choose 

¥m-m ]+t, mi-* 1)z -°- (18> 
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With four complex sources in our system, there are eight independent degrees of freedom. 
Since Z can be at most linear in each Grassmann variable, the polynomial expansion can 
be at most eighth order. All higher terms vanish due to statistics. Therefore for four sites, 
we can find an exact solution with only an eighth order Grassmann series. This polynomial 
has 69 terms which can be grouped into 14 independent symmetry classes. For example, at 
second order we have 

Z = l + ai [77(1)77(1) + 77(2)7,(2) + 77(3)77(3) + 77(4)77(4)] 
+a 2 [77(2)7/(1) + 77(3)7/(2) + 77(4)77(3) + 77(1)7/(4) 
-77(1)7/(2) - 77(2)77(3) - 77(3)77(4) - 77(4)7/(1)] 
+a 3 [77(3)7/(1) + 77(4)7/(2) + 77(1)7/(3) + 77(2)7/(4)] (19) 

where the coefficients are the Green's functions 

ai = $(1)^(1)) 

a 2 = $(2)^(1)) (20) 
a 3 = $(3)^(1)) 

Since this is an exact solution, no Galerkin procedure is needed. We can apply the Grass- 
mann equation to the full polynomial, set all the coefficients of the resulting power series to 
zero, and solve the resulting system of equations. For the exact solution, we have as many 
independent equations as unknowns and our set of algebraic equations for the coefficients 
of Z is completely specified. The result for the free two-point function is 

2M 2 + 1 



ai 



2M(M 2 + 1) 



° 2 = 2(m'+1) < 21 > 
1 

03 



2M(M 2 + 1) 

The complete results for all the Green's functions can be found in the Appendix. These 
results agree identically with these from inverting the propagator matrix, or equivalently 
from the lattice Fourier transformed solution to the Green's function equation. 
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We can generalize our action to include self-interactions 

N 1 

S = ^-^)[^ + l)-^-l)] + M^« 

i 

+ l)V(* + l) + - l)^(* - l)] 
+[^(i)+^(*M0]- (22) 

In terms of the lattice functional equations, this means adding third order derivative terms 
to the differential equations 

G S 3 Z G 5*Z 

2 5fj{i)6r}{i + l)5fj(i + 1) 2 5fj{i)6r]{i - l)5fj(i - 1) 1 ' 

We solve the interacting model just as we did for the free case. Again, we find an exact 
solution yielding the following coefficients 

2M 3 + (2G + 1)M 



ai 



2M 4 + (4G + 2)M 2 + G 2 + G 



2M 2 + G 

° 2 ~ 4M 4 + (8G + 4)M 2 + 2G 2 + 2G ( ^ 



M 

A3 



2M 4 + (4G + 2)M 2 + G 2 + G 

Notice that the coefficients respect the free field limit as g — > 0, and therefore correspond 
to the path integral solution. 

5 Galerkin Method for Fermions 

For larger lattices, even though the power series is self-terminating, it becomes too unwieldy 
for practical calculations. Therefore, we truncate the exact expansion at some lower order, 
and determine the coefficients numerically. As in the boson case, we have the problem of 
inconsistent equations. While the series expansion for Z contains the full symmetry group 
of the Action, the individual Grassmann differential operators are invariant under only 
a subgroup. The algebraic relations for the coefficients of the truncated Z are therefore 
overdetermined. This implies that an approximate procedure is needed to determine the 
coefficients. 
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The Galerkin procedure is based on an attempt to minimize the error function, or the 
residual R, about the origin in source space 



where Df is the fermionic source operator centered at site i and Zt[t], rf\ is the truncated 
partition function that depends on the Grassmann sources 77, ry. The residual represents the 
error induced by truncation. For bosons, in order to minimize the error, we defined an inner 
product in the space spanned by functions of the source variables {</(!), . . . , J(N)}. 



We demand that inner products of R with linearly independent test functions {Tj} vanish, 
and solved the resulting algebraic equations. In this way, we have fitted our approximation 
to a solution of the differential equation. Solutions obtained in the manner are called "weak 
solutions" where convergence is guaranteed in the mean. 

As a first attempt to mimic this procedure for fermions, we could begin by defining 
an inner product in the Grassmann function space 1T(J. Unfortunately, there are several 
problems with this approach. Principally, taking a Grassmann inner product of R with 
Grassmann test functions is a very severe operation. In fact, the Grassmann inner product 
kills off most of the equations that we want to solve. Furthermore, recall that for bosons, the 
Galerkin method gave average solutions in an e-neighborhood of the origin; they represented 
"weak" solutions to the differential equations. Since Grassmann integration is really a 
formal construct with no measure theoretic interpretation, then in what sense are they 
average solutions? 

Instead, we propose to sidestep all these difficulties by making a few simple observations. 
First of all, the Grassmann sources which create all of the difficulties carry virtually no 
physics. All the information about a QFT is contained in the coefficients of Z, i.e. the 
Green's functions, which are simple real numbers. The sources are arbitrary, and in fact, 
for fermionic problems, their only real role is to give the correct signs in the Green's functions 
equations for the coefficients. Once we have obtained R, all the physical information has 
been completely specified. Any manipulations of R which are mathematically consistent 
does not alter the physics. 



R = t)fZ T %r ] \ 



(25) 




(26) 
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With these observations in mind, we propose the following alternative Galerkin method 
for fermions. Once we have R, we relax the anticommutivity of all source variables, and 
allow them to be integrated in the usual way, using the boson inner product. Test functions 
are constructed as derivatives of scalar polynomials that have the same symmetry as the 
fermion invariant polynomials. The rest of the calculation proceeds identically as for the 
boson case; inner products of R with scalar test functions gives algebraic relations for the 
coefficients of Z. After solving this system of equations, e is taken to zero. 

A more compact way to view this procedure is to consider each source variable as being 
decomposed into a scalar piece times a Grassmann unit vector 



The unit vector fj(i) carries the anticommutivity and enforces the commutation relations 
while jJL acts like a length that can be integrated in the usual sense. In this framework, the 
functional equations act on the Grassmann piece while the Galerkin inner product operates 
on the scalar component. 

6 Numerical Solutions 

We test our truncation/Galerkin procedure against the four site model discussed previously. 
Since we can solve this system exactly, it is useful for checking the convergence of our 
method. Solutions for this model were found using an eighth order Grassmann polynomial. 
To examine the convergence of our method, we first truncate this series at fourth order. 
Since this is not an exact solution, we use a Galerkin approach to find an approximate one. 
The expressions for the coefficients are complicated rational functions of the parameters. 
Setting M = 1.0, the results are 



rj(i) = nfj{i). 



(27) 



2117446 G + 7294107 



ai = 



1786968 G 2 + 9742031 G + 9725476 
3092325 G + 4862738 



a 2 = 



2(1786968 G 2 + 9742031 G + 9725476) 
330478 G + 2431369 



(28) 



A3 = 



1786968 G 2 + 9742031 G + 9725476 
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g 


4th Order 


6th Order 


Exact 





0.75 


0.75 


0.75 


0.1 


0.7003 


0.7089 


0.7095 


0.5 


0.5553 


0.5840 


0.5926 


1.0 


0.4428 


0.4803 


0.5 


10 


0.0996 


0.1160 


0.1494 


mergence of fermionic Galerkin method : 


g 


4th Order 


6th Order 


Exact 





0.25 


0.25 


0.25 


0.1 


0.2413 


0.2334 


0.2328 


0.5 


0.2130 


0.1928 


0.1852 


1.0 


0.1871 


0.1677 


0.15 


10 


0.0626 


0.1141 


0.0390 



Table 2: Convergence of fermionic Galerkin method for (^(2)^(1)) 

These should be compared against the exact solution in Eq (6.78). Similarly, we truncate 
at sixth order, finding 

3 (618904 G + 510825) 



ai 
03 



1345887 G 2 + 3667541 G + 2043300 

(G + 2) (278354 G + 510825) 
2 (1345887 G 2 + 3667541 G + 2043300) 
278354 G + 510825 



(29) 



1345887 G 2 + 3667541 G + 2043300 
To see exactly how these answers converge, we substitute a range of values for G and 
examine the behavior of the coefficients. Results for the two-point functions are displayed 
in the tables. 

For large lattices, we need to perform calculations numerically. We calculated the in- 
teracting propagator for an eight site lattice in one dimension. Results are presented for a 
Grassmann polynomial truncated at fourth and sixth order. These are nontrivial calcula- 
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g 


4th Order 


6th Order 


Exact 





0.25 


0.25 


0.25 


0.1 


0.2299 


0.2223 


0.2217 


0.5 


0.1726 


0.1543 


0.1481 


1.0 


0.1299 


0.1118 


0.1 


10 


0.0201 


0.0190 


0.0065 



Table 3: Convergence of fermionic Galerkin method for (^(3)^(1)) 

tions since the fourth order polynomial alone has 124 coefficients. These coefficients must 
be determined by the Galerkin method. 

The results are checked against a mean field calculation. In general, mean field the- 
ory is not exact, but represents only a single pole approximation to the propagator. Our 
propagator, on the other hand, should include contributions from all poles. Thus, only 
loose comparisons should be made between the mean field results and the Source Galerkin 
numbers. 

It should be noted that calculations with the Source Galerkin method are very efficient, 
especially when compared to Monte Carlo. Typically, for lattice gauge theory, the fermions 
must be quenched to make calculations tractable. Here, determination of the interacting 
propagator for a system of dynamical fermions presents no special difficulty. The bulk of 
a calculation involves a single matrix inversion for a given set of parameters. This is in 
contrast to Monte Carlo where many sweeps through the lattice are necessary to reduce 
statistical error. As the plots show, the Source Galerkin calculations are very clean. They 
show rapid convergence compared to the mean field results even at intermediate couplings 
and using only low order polynomials. 



13 



4th Order 


6th Order 


Mean Field 


0.5458 


0.5831 


0.6287 


0.2288 


0.2104 


0.2218 


0.0883 


0.0753 


0.0797 


0.0177 


0.0329 


0.0243 


0.0166 


0.0267 


0.0196 



Table 4: Interacting Propagator on 8 Sites with M=1.0, g=0.5 
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Interacting Fermion Propagator 



Source Galerkin — 
Mean Field o 




Figure 1: 8 site lattice with M=1.0 and G=0.5 
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7 Leading Behavior of the Gross-Neveu Model 

Since the computational complexity grows rapidly with the number of sites, we need to 
go beyond simple polynomial expansions to examine larger systems. It is important to 
note that we are not restricted to power series expansion. Any complete set of functions is 
allowed where the choice is restricted only by convenience. We exploit this flexibility with 
respect to the choice of expansion functions. A more general choice is a Gaussian ansatz. 
This choice leaves only one unknown parameter to solve, independent of the lattice size. We 
will use this ansatz to examine the structure of a nontrivial model of interacting relativistic 
fermions. 



The Gross-Neveu model [12| is a (l+l)-dimensional renormalizable field theory of N 
species of Dirac fermions interacting through a 4-fermion term. The action is given by 

N - 1 - 

L = 53 [^ j W i] - -a 2 (V- (i ¥ (i) ) 2 ] (30) 
i=i 2 

where the four-fermion term is often replaced by a quadratic term with the help of a real 
scalar field a 

N - 1 

L = y + A <? 2 + (t^o')^i. (31) 

The action displays a discrete chiral symmetry 

V> -> -^75 
a — > —a 

where 75 = «7i72- The model can be studied analytically by the saddle point method, 
leading to a 1/N expansion for g 2 N fixed. At leading order, the chiral condensate is given 
by 

(<j)=Ae-^ 2 ^. (32) 

The theory describes a free fermion where the dynamical mass nif = (a) results from 
spontaneous chiral symmetry breaking. 
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We apply the Source Galerkin method to this model. For simplicity, we employ stag- 
gered fermions. To probe the details of the system, we need lattices considerably larger 
than the ones used in the ID example. For such lattices, construction and manipulation 
of the Grassmann invariant polynomials becomes difficult. Instead, we take advantage of 
the extreme flexibility of the Galerkin method with regards to the choice of expansion func- 
tions. Following a similar calculation for scalar ||], we determine the leading asymptotic 
behavior for vanishing sources J — > of the GNM. 

After including sources (fj n ,r) n , S n ) for the fermion and auxiliary fields, we vary the 
staggered GN action to obtain the equations of motion 



(£> n , m + <7„) Xn = Vn 



g + E*nXn = S n (33) 
n=l 



where 



D n , m = ~ [S U;m+ i ~ S nim _i\ + ^ (-) ni [5„ )m+ 2 - 5„ im _2] C 34 ) 

is the Susskind differential operator with 1,2 being unit vectors along the two Euclidean 
directions and where A = g 2 N . After taking matrix elements and inserting functional 
definitions, we obtain the lattice functional equations 

1 SZ « 5 2 Z 
We solve these equations with a Gaussian in the sources 

% — e W(fj,r]) 

and find the following equations for W{f) } rf) 



(36) 



SW6Z_ 
1 n ' m+ 5S n l 5fj- ~ Vn 

N 5W ^ r 5 2 W 5W \,SW 

Expanding W in a 2nd order power series in the Grassmann sources gives 

W(fj,r!) = fi k G(k-j) Vj (38) 
17 



where G(k — j) is the interacting propagator. Since a n is not dynamical, we ignore it at 
this order, and set SW/5S which measures chiral symmetry breaking equal to a constant 

SW 

— = (a) = a . (39) 
Keeping only linear terms in the sources, we obtain a dynamical equation 

[D n ,m + cr ] G(n - m) = 5 n , m (40) 
as well as a constraint equation. 

j = G(0). (41) 

The dynamical equation suggests that G(n — m) has the form of a free staggered propagator 
with mass ctq. We can solve this equation by Fourier transform 

G(n - m) = -L V V ^o + ^in(pi) + ,(-)^ sin(p 2 ) &ip{n _ m) 
L " ^ 0p ^ ^+shi 2 (pi)+sin 2 (p 2 ) 

where n, m,p are two dimensional vectors. Putting this expression into the constraint gives 

a gap equation for gq that must be solved self-consistently 

1 1 1 

~ X = T2 pSpS -o 2 + sm 2 (pi)+sin 2 (p 2 )- (43) 
This type of equation is familiar from large ./V expansions. 

It is important to note that this solution has been obtained by fitting a Gaussian 

Z = expw(f?, rj) (44) 

to the functional equations where we kept only second order terms in the expansion for W 

w(v,v) = VjGj,kVk + ■■■ (45) 

At this level of approximation, the procedure may seem rather trivial. The results we have 
obtained are closely related to the usual large N expansion results. However, it should be 
noted that our numerical procedure has determined the leading order solution of a nontrivial 
fermionic problem. Systematic improvements of approximations are then possible where the 
Gaussian is only the leading contribution. By considering Z as a Gaussian multiplied by a 
power series, it is possible to improve our leading order result. The Galerkin method can 
be used to determine the coefficients of the polynomials. 

18 



8 Conclusion 



We have presented a new numerical method to solve lattice fermion theories. It is based 
on the differential formulation of QFT in the presence of an external source. By examining 
the functional differential equations for a theory on a finite lattice, we obtained a set of 
coupled Grassmann PDEs for the generating functional Z. For nonlinear field theories with 
polynomial interactions, these equations are always linear. Once Z is obtained, the lattice 
Green's functions can be extracted by functional differentiation. 

To construct solutions, we can expand Z on any complete set of functions in the source 
variables {rji}. A particularly simple choice is polynomial functions. We saw that these 
functions gave very rapid convergence even using low order polynomials. Calculations were 
efficient to perform and produced very clean numbers. The bulk of any calculation involved 
only a single matrix inversion. The matrix to be inverted becomes large, but unlike Monte 
Carlo, there is only one inversion for a given set of parameters. 

Due to computational complexity, polynomial basis functions are limited to small lat- 
tices. For larger systems, more general schemes are possible. The Galerkin procedure 
of fitting approximate solutions to the functional equations allows tremendous freedom in 
choosing expansion functions. We illustrated one possibility by using a Gaussian to extract 
the leading behavior of the solution. In this formulation, we had only one unknown param- 
eter, independent of the size of the lattice. Systematic corrections can then be calculated. 
We will have more to say about other generalized approaches in future communications. 

This method is deterministic, and therefore, does not suffer from many of the problems 
usually associated with numerical solutions of fermionic problems. For example, Monte 
Carlo simulations of lattice gauge theories have difficulty including dynamical fermions. 
But with the Source Galerkin method, calculation of the interacting propagator presents 
no special complication. In addition, fermionic systems at nonzero chemical potential have 
problems defining a positive definite probability measure for simulations due to the "minus 
sign" problem. Since our method is deterministic, this is not an issue. 
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9 Appendix A: Exact Coefficients for Four Site Model 



2 M 2 + (2 G + 1) M 



A2 = 



2 M 4 + (4 G + 2) M 2 + G 2 + G 
2M 2 G 

2M 4 + (4G + 2)M 2 + G 2 + G 



M 

A3 = 



A4 = 



2M 4 + (4G + 2)M 2 + G 2 + G 

4M 2 + 2G + 1 
4M 4 + (8G + 4)M 2 + 2G 2 + 2G 



A5 M 



A6 



2M 4 + (4G + 2)M 2 + G 2 + G 



4 M 4 + (8 G + 4) M 2 + 2 G 2 + 2 G 



M 

A7 = 



A8 = 



A9 = 



2M 4 + (4G + 2) M 2 + G 2 + G 
1 

4 M 4 + (8 G + 4) M 2 + 2 G 2 + 2 G 
2M 2 

2M 4 + (4G + 2)M 2 + G 2 + G 



A10 = 



All 2M 



A12 



2M 4 + (4G + 2) M 2 + G 2 + G 



2M 4 + (4G + 2) M 2 + G 2 + G 



A13 = 

20 
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Appendix B: Fourth Order Coefficients with M 

2117446 G + 7294107 



Al 



A2 = 



A3 = 



A4 = 



A5 



1786968 G 2 + 9742031 G + 9725476 

3092325 G + 4862738 
2(1786968 G 2 + 9742031 G + 9725476) 

330478 G + 2431369 
1786968 G 2 + 9742031 G + 9725476 

3573936 G + 12156845 
2(1786968 G 2 + 9742031 G + 9725476) 

792533 G + 2431369 
1786968 G 2 + 9742031 G + 9725476 



2431369 

A6 ■■ 



A7 = 



A8 = - 



A9 = 



2(1786968 G 2 + 9742031 G + 9725476) 

792533 G + 2431369 
1786968 G 2 + 9742031 G + 9725476 

736530 G- 2431369 



2(1786968 G 2 + 9742031 G + 9725476) 
1061159 G + 4862738 



AW 



1786968 G 2 + 9742031 G + 9725476 
24642 G 

1786968 G 2 + 9742031 G + 9725476 
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Appendix C: Sixth Order Coefficients with M 



3 (618904 G + 510825) 



A2 



A3 = 



A4 



A5 



A6 



A7 



A8 



A9 = 



1345887 G 2 + 3667541 G + 2043300 

(G + 2)(278354 G + 510825) 
2(1345887 G 2 + 3667541 G + 2043300) 

278354 G + 510825 

1345887 G 2 + 3667541 G + 2043300 

5(482684 G + 510825) 
2(1345887 G 2 + 3667541 G + 2043300) 

278354 G + 510825 
1345887 G 2 + 3667541 G + 2043300 

278354 G + 510825 
2(1345887 G 2 + 3667541 G + 2043300) 

278354 G + 510825 
1345887 G 2 + 3667541 G + 2043300 

278354 G + 510825 
2(1345887 G 2 + 3667541 G + 2043300) 

2(278354 G + 510825) 
1345887 G 2 + 3667541 G + 2043300 



A10 = 



1021650 

All ■- 



A12 = 



A13 = 



1345887 G 2 + 3667541 G + 2043300 
510825) 

1345887 G 2 + 3667541 G + 2043300 

556708 G 

1345887 G 2 + 3667541 G + 2043300 
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